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ABSTRACT 

Using Herschel data from the deepest SPIRE and PACS surveys (HerMES and 
PEP) in COSMOS, GOODS-S and GOODS-N, we examine the dust properties of in- 
frared (IR)-luminous (Lir>10^° Lq) galaxies at 0.1 < z < 2 and determine how these 
evolve with cosmic time. The unique angle of this work is the rigorous analysis of 
survey selection effects, making this the first study of the star-formation-dominated, 
IR-luminous population within a framework almost entirely free of selection biases. We 
find that IR-luminous galaxies have spectral energy distributions (SEDs) with broad 
far-IR peaks characterised by cool/extended dust emission and average dust temper- 
atures in the 25-45 K range. Hot (T > 45) SEDs and cold (T < 25), cirrus-dominated 
SEDs are rare, with most sources being within the range occupied by warm starbursts 
such as M82 and cool spirals such as M51. We observe a luminosity-temperature {L—T) 
relation, where the average dust temperature of log [Lir/L0]=12.5 galaxies is about 
10 K higher than that of their log [Lir/L0]=1O.5 counterparts. However, although the 
increased dust heating in more luminous systems is the driving factor behind the 
L — T relation, the increase in dust mass and/or starburst size with luminosity plays a 
dominant role in shaping it. Our results show that the dust conditions in IR-luminous 
sources evolve with cosmic time: at high redshift, dust temperatures are on average up 
to 10 K lower than what is measured locally (2 < 0.1). This is manifested as a flatten- 
ing of the L~T relation, suggesting that (U)LIRGs in the early Universe are typically 
characterised by a more extended dust distribution and/or higher dust masses than 
local equivalent sources. Interestingly, the evolution in dust temperature is luminos- 
ity dependent, with the fraction of LIRGs with T <35K showing a 2-fold increase 
from z ■^O to z ~2, whereas that of ULIRGs with T <35K shows a 6-fold increase. 
Our results suggest a greater diversity in the IR-luminous population at high redshift, 
particularly for ULIRGs. 



1 INTRODUCTION 

The discovery of a class of infrared (IR)-luminous 
(LiR > 10'^'' Lq) galaxies in the 60s (e.g. Johnson 1966; Low 
& Tucker 1968; Kleinmann & Low 1970), followed by the de- 
tection of the cosmic infrared background (Puget et al. 1998; 
Hauser et al. 1998) unfolded a new era in extragalactic as- 
tronomy. Infrared/submm surveys with IRAS (Neugebauer 
et al. 1984), ISO (Kessler et al. 1996), JCMT/SCUBA (Hol- 
land et al. 1999), Spitzer (Werner et al. 2004) and AKARI 
(Murakami et al. 2007), revealed that the early Universe 
was more active than previously thought, uncovering a large 
number of dust-enshrouded galaxies whose bolometric lumi- 
nosity emerges almost entirely in the infrared (e.g. Soifer et 
al. 1984a; Sanders & Mirabel 1996; Lutz et al. 1996; 2005; 
Rowan-Robinson et al. 1997; 2005; Lisenfeld, Isaak & Hills 
2000; Goto et al. 2010 and many more). These IR-luminous 
galaxies are rare in the local Universe (e.g. Kim & Sanders 
1998) but exhibit a strong increase in number density at ear- 
lier epochs (e.g. Takeuchi, Buat & Burgarella 2005), being 
responsible for about half the total light emitted from all 
galaxies integrated through cosmic time (e.g. Gispert et al. 
2000; Lagache et al. 2005; Dole et al. 2006). The abundance 
of these sources at high redshifts {z ~l-4; e.g. Hughes et al. 
1998; Bales et al. 1999, 2000; Blain et al. 2004; Le Floc'h et 
al. 2004; Schinnerer et al. 2008; Pannella et al. 2009) indi- 
cates that they started forming very early in cosmic history, 
potentially challenging the hierarchical paradigm of ACDM 
(e.g. Granato et al. 2004; Baugh et al. 2005; Bower et al. 
2006; Somerville et al. 2008). 

IR-luminous galaxies are the ideal laboratories for stud- 
ies of galaxy formation and evolution through chemical en- 
richment, star-formation, black hole accretion and stellar 
mass build-up. They hide an immensely active interstel- 



lar medium (ISM; e.g. Lutz et al. 1998; Farrah et al. 2003; 
Narayanan et al. 2005; Sturm et al. 2010) and are the ulti- 
mate stellar nurseries, with star-formation rates (SFRs) up 
to a few thousand times higher than Milky Way (MW)-type 
galaxies (e.g. Kennicutt 1998; Egami et al. 2004; Choi et al. 
2006; Rieke et al. 2009). In addition, they are amongst the 
most massive galaxies in the Universe (e.g. Dye et al. 2008; 
Michalowski et al. 2010) and often their morphologies show 
signs of interactions and mergers (e.g. Sanders & Mirabel 
1996; Farrah et al. 2001; 2002; Moustakas et al. 2004; Kar- 
taltepe et al. 2010b). Finally, they frequently harbour an 
active galactic nucleus (AGN), which is commonly consid- 
ered a key player in the evolution of the system (e.g. Genzel 
et al. 1998; Ptak et al. 2003; Alexander et al. 2005; Page et 
al. 2012). 

Until recently our view of the IR-luminous galaxy 
population at high redshift has been based on data 
from the space observatories 75*0 and Spitzer, as well as 
ground-based submm/mm facilities such as JCMT/SCUBA, 
APEX/LABOCA, IRAM/MAMBO and SMA/AzTEC. Ah 
though huge advances have been made in our understanding 
of the nature and evolution of these sources, it has been chal- 
lenging to reconcile the data from space observatories with 
comparable ground-based IR/mm datasets. In recent years 
it has become increasingly apparent that, besides strong evo- 
lution in IR galaxy number density (e.g. Le Floc'h et al. 
2005; Huynh et al. 2007; MagneUi et al. 2009, Berta et al. 
2010; 2011), the physical properties of IR galaxies might 
also evolve with time, with the rate of evolution potentially 
changing as a function of luminosity (Seymour et al. 2010). 
Studies of the local Universe showed that ultraluminous in- 
frared galaxies (ULIRGs) are characterised by warm aver- 
age dust temperatures (e.g. Soifer et al. 1984; Klaas et al. 
1997; Clements, Dunne & Eales 2010), strong silicate ab- 
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sorption and PAH emission features in their mid-IR contin- 
uum (e.g Brand! et al. 2006; Armus et al. 2007), as well as 
compact starburst sizes (e.g. Condon et al. 1991; Soifer et 
al. 2001). However, with the onset of submm/mm facilities 
which probed the early Universe {z >2), such as SCUBA 
in the late 1990s, a different picture emerged. Many IR- 
luminous galajdes at high redshift were found to be less com- 
pact than their local counterparts (e.g. Tacconi et al. 2006; 
lono et al. 2009; Rujopakarn et al. 2011), exhibiting stronger 
PAH emission (Farrah et al. 2007; 2008; Vahante et al. 2007) 
and a greater abundance of cold dust (Kovacs et al. 2006; 
Pope et al. 2006; Coppin et al. 2008). It was later shown 
that these differences in the measured dust properties were 
partly due to selection effects and partly due to evolution 
(Symeonidis et al. 2009; 2011a). Moreover, the exploitation 
of long wavelength data from ISO (Rowan- Robinson et al. 
2005) and Spttzer (Symeonidis et al. 2007; 2008), enabled 
the discovery of IR-luminous galaxies at <1, with a spec- 
trum of properties which overlapped with both the local 
population and the SCUBA-detected, z >2, sources, pro- 
viding the missing link between the two (Symeonidis et al. 

2009) . 

The launch of the Herschel Space Ohservator-^ (Pil- 
bratt et al. 2010) has opened a new window in infrared 
astronomy, as it is the only facility to date and for the 
foreseeable future to perfectly span the wavelength range 
in which most of the Universe's obscured radiation emerges 
(70-500 /xm), uncovering unprecedented numbers of dust- 
enshrouded galaxies over a sizeable fraction of cosmic time. 
The large dynamical range of the PACS (Poglitsch et al. 

2010) and SPIRE (Griffin et al. 2010) instruments, have en- 
abled spectral energy distributions (SEDs) to be compiled 
for a large range of objects, both for AGN (e.g. Hatzimi- 
naoglou et al. 2010; Seymour et al. 2011) and star-forming 
galaxies (e.g. Rowan-Robinson et al. 2010). Recent studies of 
the properties of the IR-luminous galaxy population using 
Herschel data provide an excellent showcase of the capa- 
bilities of this observatory (some examples from the multi- 
tude of Herschel papers on this topic: Amblard et al. 2010; 
Rowan-Robinson et al. 2010; MagneUi et al. 2010; 2012, 
Magdis et al. 2010; 2011; Gruppioni et al. 2010; Eales et al. 
2010b; Hwang et al. 2010; Elbaz et al. 2010; 2011; Rodighiero 
et al. 2010; Ivison et al. 2010; Buat et al. 2010; Dye et al. 
2010; Berta et al. 2010; 2011; Dunne et al. 2011; Symeoni- 
dis et al. 2011b; Kartaltepe et al. 2012). Results from these 
studies carried out during the Science Demonstration Phase 
(SDP) of the largest extragalactic surveys, HerMES [Her- 
schel multi-tiered extragalactic survey; Oliver et al. 2012), 
PEP (PACS Evolutionary Probe; Lutz et al. 2011) and H- 
ATLAS {Herschel Astrophysical Terahertz Large Area Sur- 
vey; Eales et al. 2010a), confirmed previous findings on the 
diversity of IR-luminous galaxy properties at high redshift, 
as well as the existence of high-redshift sources with no local 
equivalents, moving us closer in understanding the complex 
nature of the infrared galaxy population up to z ~3. 

In this paper we report a comprehensive study of the 
SEDs and dust temperatures of IR-luminous galaxies up 



^ Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and 
with important participation from NASA. 



to z ^2, using the deepest available Herschel /VACS and 
SPIRE data acquired as part of the PEP and HerMES con- 
sortia in the COSMOS and GOODS (N & S) fields. A key 
aspect of our work is the attempt to eradicate survey selec- 
tion effects, an issue which has plagued previous attempts 
to canvas the range of infrared SEDs (see Symeonidis et al. 
2011a). Thus for the first time we are able to examine the 
aggregate properties (infrared luminosity, dust temperature, 
SED shape) of IR galaxies within an almost entirely bias-free 
framework. The paper is laid out as follows: the introduc- 
tion is followed by a section on the sample selection (section 

01 and SED measurements (section [3} . In section 3] we dis- 
cuss how we deal with AGN, in order to obtain a sample 
which is star-formation dominated in the infrared. Section 
Ois devoted to treatment of selection effects, enabling us to 
assemble a complete sample of IR galaxies. In section |6] we 
present our results and discuss them in section [T] Finally 
our summary and conclusions are presented in section [S] 
Throughout we employ a concordance ACDM cosmology of 
Ho=70kms"^Mpc"^ r2M=l-r2A=0.3. 

2 THE HERSCHEL SAMPLE 
2.1 Initial selection 

The starting point for this work are data from Herschel, cov- 
ering three extragalactic fields: the Great Observatories Ori- 
gins Deep Survey (GOODS)-North and South (Giavalisco 
et al. 2002) and the Cosmic Evolution Survey (COSMOS) 
field (Scoville et al. 2007). We use PACS 100 and 160 /xm 
and SPIRE 250, 350 and 500 /^m images, acquired as part of 
PEP and HerMES respectively. Source extraction from the 
PACS and SPIRE imagef0 is performed on the IRAC-3.6 ^.m 
positions of the f24^30/iJy GOODS (N and S) sources and 
f24^60^Jy COSMOS sources, as described in MagneUi et 
al. (2009) and Roseboom et al. (2010; 2012). This method 
of source extraction on prior positions is widely used and 
enables identifications of secure counterparts over the whole 
SED. In this case however, its significant advantage, lies in 
its ability to effectively deal with source blending in the 
Herschel bands, particularly for SPIRE where the beam is 
large (18.1, 24.9 and 36.6 arcsec FWHM at 250, 350 and 
500 /im respectively; Nguyen et al. 2010). By using prior in- 
formation to identify galaxies in the Herschel images, we are 
able to extract 'clean' photometry for each galaxy, even for 
those which appear blended in the PACS and SPIRE bands. 
For information on the GOODS Spiteer/MIPS 24 /^m data 
see MagneUi et al. (2009) and for information on the COS- 
MOS Spitzer/Ml¥S 24 ^im data see Sanders et al. (2007); Le 
Floc'h et al. (2009). The 3cr sensitivity limits of the PACS 
100 and 160 /im catalogues respectively are 5 and 10 mjy 
for COSMOS, 3 and 6mJy for GOODS-N and 1 and 2mJy 
for GOODS-S. A 3 a detection in SPIRE using prior posi- 
tions and the cross-identification method of Roseboom et al. 
(2010) is approximately 8, 11 and 13mJy at 250, 350 and 
500 /im in all fields. In the case of the PACS bands a is only 
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Table 1. Table showing the detection statistics of the Herschel sample used in this work. The first and second columns show the field 
and number of 24 fim sources with /24>30 /ijy and /24>60 fijj for GOODS (N & S) and COSMOS respectively, whose positions are used 
as priors for source extraction in the Herschel bands. Column 3 shows the total number of 'isolated' 24 ^m sources defined as having no 
companion within 8"; the percentage in brackets is calculated out of the number in column 2. Column 4 shows the fraction detected at 
100 and 160 ^m, whereas column 5 shows the fraction detected at 160 and 250 /xm. The final column shows what fraction of the 24 /xm 
population is detected when the two criteria are used in disjunction, i.e. [100 and 160 /xm] OR [160 and 250 fim]. This is the criterion 
used to select the initial Herschel sample f section 12.111 . The fractions shown in columns 4, 5 and 6 are out of the number of sources 
in column 3. As also mentioned in section [2.11 the Scr sensitivity limits of the PACS 100 and 160 /xm catalogues respectively are 5 and 
lOmJy for COSMOS, 3 and 6mjy for GOODS-N and 1 and 2mjy for GOODS-S. A 3cr detection in SPIRE using prior positions and 
the cross-identification method of Roseboom et al. (2010) is approximately 8, 11 and 13mjy at 250, 350 and 500 /im in all fields. 



field 


number of 


number of 


fraction detected at 


fraction detected at 


fraction detected at 




24/im sources 


'isolated' 24 /im sources 


100-fl60Aim (>3cr) 


160-1-250 Aim {>3a) 


[100+160 Aim] OR [160+250 Mm] 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


GOODS-N 


2149 


1401 (65%) 


7% 


7% 


9% 


GOODS-S 


2252 


1580 (70%) 


21% 


9% 


22% 


COSMOS 


52092 


33407 (64%) 


5% 


6% 


7% 



the photometric error, whereas for the SPIRE bands, a in- 
cludes confusion error (see Nguyen et al. 2010 for the SPIRE 
confusion limits). 

The initial selection for our sample includes all 24 fj.m 
sources that have detections (at least 3a) at [100 and 
160 Aim] OR [160 and 250 Aim] (where 'OR' is the oper- 
ator representing disjunction in Boolean logic). This en- 
sures that (i) the sample consists of IR-luminous galaxies 
(LiR > 10^" Lq), (ii) the sample is as complete as possible 
over a large redshift range with respect to SED types, given 
the PACS and SPIRE selection functions (see section [5] for 
more details) and (iii) there are at least 3 reliable photo- 
metric points in the SED (24Atm+ 2 Herschel bands) for 
subsequent measurements. 

Our selection is in essence a colour selection rather than 
a single band selection, as we require sources to be detected 
at both 24 and 160 Aim (the additional Herschel photometry 
at 100 or 250 Atm has a small effect on the sample selection 
but enables more accurate SED measurements; see section 
[5l. As a result, the detection rate is more strongly dependent 
on the SED shape, in our case the mid-to-far-IR continuum 
slope, than for a typical flux limited survey. Given the flux 
limits reported earlier, the 24 Atm survey is 66, 166 and 200 
times deeper than the PACS 160 Atm survey in GOODS-S, 
COSMOS and GOODS-N respectively. We find that the dif- 
ferent ratios in flux density limits (/J6o//24") between the 3 
fields, introduce a bias with respect to the SED shapes that 
are detected in each survey particularly for objects with flux 
densities close to the limit. To mitigate this effect, we match 
the relative PACS-24 Aim depths of the GOODS flelds to the 
COSMOS survey, as the latter covers the largest area and 
hence dominates the statistics of the final sample. This is 
done as follows: the GOODS-S sample is cut at /i6o=5mJy 
and the GOODS-N sample is cut at /24 =36AiJy, so that 
/{go 7/24" ^ 166 in all cases. Note that matching the sam- 
ples in this way ensures that the relative biases between the 
surveys are minimised, i.e. that all three surveys probe the 
same range of SED types, however it does not deal with 
absolute biases; these are dealt with in section [S] 

Besides the photometric selection criteria, we also re- 
strict the sample to sources which have no other 24 Aim com- 
panions within 8"; i.e. 'isolated' sources. This allows us to 
work with more reliable photometry, as at the longer wave- 



lengths, where the Herschel beam is large, flux extraction in 
the Herschel bands can be problematic when dealing with 
blended sources. The choice of an 8"radius is larger than the 
24 Atm beam (6") and visual inspection shows that it is suf- 
ficient to eliminate problematic blends. In addition, 8"is the 
scale of the first airy ring of bright 24 Atm sources. In cases 
where there is a companion within the first airy ring, the 
Herschel flux is assigned to the bright source in the centre, 
if the companion source is 10 or more times fainter (Rose- 
boom et al. 2010). By eliminating such cases from our sam- 
ple, we avoid the occurance of a potential bias causing the 
measured Herschel flux density to positively correlate with 
the 24 Aim flux density. As a result, for the 24 Aim sources 
we subsequently use, there have been no prior assumptions 
when assigning Herschel fluxes. This is confirmed by per- 
forming a K-S test on the flux density distribution of the 
whole 24 Aim population and that of the 'isolated' 24 Atm 
sources in each field. We find the two to be entirely consis- 
tent, suggesting that our approach works well in eliminating 
problematic sources without introducing systematic biases. 

At this point, our assembled Herschel sample consists 
of 2500 sources (2206 from COSMOS, 173 from GOODS- 
S and 121 from GOODS-N). Some statistics for the initial 
sample are shown in table [T] 



2.2 Redshifts 

The redshifts we use are a combination of spectroscopic and 
photometric, assembled from various catalogues: Berta et 
al. (2011) for GOODS-N, Cardamone et al. (2010) and San- 
tini et al. (2009) for GOODS-S and Ilbert et al. (2009) and 
Salvato et al. (2009); (2011) for COSMOS. The optical po- 
sitions of sources in these catalogues are cross-matched to 
the 24 Atm positions within 1". The excellent photometric 
coverage of these fields and high quality photometric red- 
shifts available, result in >90 per cent of the sources in our 
sample having a usable redshift (25 per cent spectroscopic; 
although for the final sample about 1/3 are spectroscopic), 
leaving 2313 Herschel sources for subsequent analysis. For 
more details on the quality of photometric redshifts see Ap- 
pendix]^ where we also show that the use of photometric 
redshifts does not bias our results. 
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Figure 1. Total infrared luminosity as a function of redshift for 
the initial Herschel sample of 2313 sources used for SED mea- 
surements (section [Sj . 

3 MEASUREMENTS 

We fit the photometry of the Herschel sample, with the 
Siebenmorgen & Kriigel (2007, hereafter SK07) library of 
7208 models, built on the formulation of Kriigel & Sieben- 
morgen (1994). This library of templates ranges in 5 free 
parameters, in physically acceptable combinations and the 
shape of each template is determined by the combination of 
parameters which define it: 

• the radius of the IR emitting starburst region [R) , tak- 
ing discrete values of 0.35, 1, 3, 9, 15 kpc 

• the total luminosity of the system (Ltot) ranging from 
lO^'^ to 10^^ Lq 

• the visual extinction from the edge to the centre of the 
starburst, taking discrete values of 2.2, 4.5, 7, 9, 18, 35, 70 
and 120 mag 

• the ratio of the luminosity of OB stars with hot spots 
to the total luminosity {Los/Ltot) taking discrete values of 
40, 60 and 90 per cent for the ^C3 kpc models and 100 per 
cent for the 9 and 15kpc models 

• the dust density in the hot spots in units of hydrogen 
number density (cm~^) ranging from 100 to 10000, in dis- 
crete steps 

As mentioned earlier, 6-band photometry is available — 
24 ^im from Spiteer/MIPS, 100, 160 from Herschel/ PACS 
and 250, 350, 500/im from Herschel /SPIRE — and at least 
3 bands, always including the 24 and 160 jj,m data, are used 
in the fitting. The normalisation of the templates is varied in 
order to minimise ■ The 0.68 lower and upper confidence 
limits for our computed parameters resulting from the fits 
(e.g. total infrared luminosity etc.) are calculated according 
to the prescribed confidence intervals for one interesting 
parameter, namely Xmin + li where Xmin is the minimum x^- 

Note that because our photometry only sparsely sam- 
ples the full IR SED, the parameters that characterise the 



best fit SEDs within Xmin + 1 are often degenerate and not 
well constrained. In addition, each object in the sample is 
fit with the entire SK07 library irrespective of the inherent 
luminosity of the templates; the total infrared luminosity, 
LiR (Lq), of each source is computed by integrating the 
best matched SED model between 8 and lOOO/xm (see Fig. 
[T]for a plot of LiR as a function of redshift). This allows us 
to stay clear of any assumptions which link the SED shape 
to the luminosity. However, it also implies that we cannot 
directly use some of the SK07 parameters to describe our 
sample, as they require scaling. One such parameter is the 
starburst size, as its impact on the SED shape depends on 
the total input luminosity. Finally, the SK07 grid, although 
more flexible than most stand-alone SED libraries currently 
in the public domain, is still too coarse to allow complete 
characterisation of the physical properties of the sample. 
For these reasons, we opt to use one parameter to describe 
the overall shape of the SK07 SED templates; we refer to 
this as the flux (J-"), calculated as log [Ltot /47r_R^] in units 
of Lq kpc~^, where Ltot is the given luminosity of the tem- 
plate (not our computed Lir) and R is the starburst size 
that corresponds to that template. In the SK07 formula- 
tion, for constant ^v, as R becomes larger, the dust mass 
increases as a function of R^ and hence becomes cooler, with 
the temperature being a function of Ltot/R?- Hence, the 
larger F is, the more flux reaches the edge of the starburst 
region and therefore the dust emission is warmer. One can 
interpret high and low values of F as representative of sys- 
tems with warmer/more compact and cooler/more extended 
dust-emission respectively. 

In order to calculate dust temperatures, we use a 
modified black-body function (a grey-body), of the form 
B\{T){l — e^'^^), with a wavelength dependent optical depth 
{lm^J.m/\Y (e.g. see Klaas et al. 2001) and a 
dust emissivity index /3. We assume a low opacity limit, so 
approximate the term (1 — e~^-^) by A"''. Typical reported 
values of /3 range between 1.5-2 (e.g. Dunne et al. 2000, 
Lisenfeld, Isaak & Hills 2000) and we adopt /3=1.5, con- 
sistent with studies of the far-IR emissivity of large grains 
(Desert, Boulanger & Puget 1990). The temperatures are 
derived by fitting all photometry at A^Ai„j,^_i, where i 
denotes a Herschel band (100, 160, 250, 350, 500 /xm) and 
imax is the band which corresponds to the maximum flux 
{vfv)', the 24 /im photometry is never included in the fit- 
ting. Note that although some studies have shown that a 
two-temperature greybody model (e.g. Klaas et al. 2001; 
Dunne & Eales 2001), is a more accurate description of far- 
IR dust-emission, this would not work with our available 
photometry; much longer wavelength data would be needed 
especially for sources at high redshift. However, in any case, 
our aim is to measure the average dust temperature of the 
far-IR peak for each source, thus a single temperature grey- 
body model is required; our method gives a temperature 
which is most representative of the peak dust emission. 



4 DEALING WITH AGN CONTAMINATION 

As this work targets the properties of the star-forming 
galaxy population, objects whose infrared energy budget po- 
tentially includes a significant contribution from an AGN 
need to be removed from the sample. Although the frac- 
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Figure 2. Observed Spitzer /TR,AC colours {fs.s/fs.e versus 
fs/f4.5 ) for the Herschel sample (black crosses). Red circles 
denote sources which are identified as AGN-dominated in the 
near/mid-IR using the IRAC criteria outlined in Donley et al. 
(2012). 



tion of Herschel galaxies found to host AGN is high (~30 
per cent; Symeonidis et al. 2011b), it has been shown that 
AGN in far-IR-selected galaxies do not often dominate the 
infrared or total energy budget of the system. According 
to an energy balance argument, if the AGN is energetic 
enough to contribute significantly to the infrared emission 
of a starburst galaxy, then its signature is likely to emerge 
in the mid-IR part of the SED in the form of a power-law 
continuum (e.g. Symeonidis et al. 2010). As a result, the 
most suitable way to identify such objects is by examin- 
ing their colours in the Spitzer /IRAC (3.6, 4.5, 5.8, 8 ^m) 
bands. Until recently, the most commonly used IRAC AGN 
selection criteria have been those presented in Lacy et al. 
(2004) and Stern et al. (2005). However, as shown in Yun et 
al. (2008) and Donley et al. (2012), there is a non-negligible 
chance that IR/submm selected galaxies will be erroneously 
identified as AGN-dominated in the IRAC bands. In addi- 
tion, Hatziminaoglou et al. (2009) reported 'cross-talk' be- 
tween the AGN and SB loci in the Lacy et al. (2004) dia- 
gram. Here, we use the Donley et al. (2012) IRAC criteria, 
shown to be effective in picking out AGN and sufficiently ro- 
bust against misidentifications. Fig.[2]shows /s.s/Ja.e colour 
against fs/fi.s colour for the Herschel sample. Indicated in 
red are the sources which satisfy the Donley et al. (2012) 
criteria and are hence classified as AGN-dominated in the 
near/mid-IR. These 87 (out of 2313) objects, '-^4 per cent, 
are subsequently excluded from the sample, leaving 2226 
sources. Note that this is not the fraction of AGN hosted by 
Herschel sources, rather it is the fraction of objects where 
the AGN could contribute substantially in the mid-infrared 
and hence interfere with our analysis. A 4 per cent fraction 
of sources with AGN-dominated near/mid-IR SEDs is in line 
with results from Symeonidis et al. (2010) and Pozzi et al. 



(2012) who show that AGN rarely contribute more than 20 
per cent in the IR emission of far-IR detected systems (see 
also Hatziminaoglou et al. 2010; Page et al. 2012; Rosario et 
al. 2012; Nordon et al. 2012). 



5 DEALING WITH SELECTION EFFECTS 
5.1 Herschel selection 

To accurately characterise the aggregate properties of the 
IR- luminous population and their evolution with redshift, 
there should be no bias with regard to the SED types we can 
observe, particularly regarding the far-IR where all our mea- 
surements are performed. As a result, the accuracy of our 
work rests on minimising selection biases and assembling a 
sample within an unbiased part of the L — z parameter space. 
For this purpose, we use the method described in Syme- 
onidis et al. (2011a) to examine the selection functions of 
the PACS and SPIRE bands, mapping out an SED-redshift- 
luminosity parameter space at the fiux density limits of the 
GOODS and COSMOS surveys used in this work. 

Fig. [3] shows the Lir^IO^^L© and 10^^ L0 selection 
functions for MIPS 24 fim, PACS 100, 160 ^m and SPIRE 
250, 350 and 500 fxm at the COSMOS fiux density limits. 
As also explained in detail in Symeonidis et al. (2011a), Fig. 
[3] is created by using all SED templates from the SK07 li- 
brary, normalising them to the required total IR luminos- 
ity and then scaling and redshifting them to the observed 
frame. Each SED template is then convolved with the MIPS, 
PACS and SPIRE filter transmission curves in order to get 
the weighted integrated fiux within each filter. We subse- 
quently perform a colour correction (according to the pre- 
scription in the instruments' observer manuals) in order to 
obtain a monochromatic fiux density in each band, derived 
with the same spectral shape used to calculate the measured 
flux density of real sources. For each band we then compare 
our template monochromatic flux density to the flux den- 
sity limit, in order to determine whether an object with the 
given redshift, luminosity and SED shape would be part of 
our sample. This results in the selection functions presented 
in Fig. [31 Red thick patterns mark the regions where all 
templates of a given peak wavelength are detected, whereas 
the blue and green dashed patterns indicate regions where 
only 90 and 70 per cent of the SK07 templates are recovered 
respectively. The detection rate relates to variations in SED 
shape; for example, for the MIPS 24 /^m selection function, 
only luminous SEDs with strong PAH features will be de- 
tected at z ^^1.7. In the unshaded areas less than 70 per cent 
of templates are detectable, a fraction which reduces to zero 
above a certain redshift. 

Note that the selection functions of PACS and SPIRE 
overlap significantly, suggesting that the SEDs of most 
sources will be fully sampled by Herschel photometry. How- 
ever, there are large differences between the MIPS/24 /xm 
and the Herschel selection functions, especially at 500 fim. 
This is not surprising as they probe different parts of the 
SED and it is expected that the long wavelength one even- 
tually turns to favour very cold sources and the short wave- 
length one eventually turns to favour warm sources. The 
24^m selection is a steep function of SED shape: SEDs with 
a significant warm dust component are favoured up to very 
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Figure 3. log (Lir/L0)=11 and log (Lih/L0)=12 selection func- 
tions for PACS and SPIRE, constructed using the SK07 model 
library; see also Symeonidis et al. (2011a). The plot shows SED 
peak wavelength (left y-axis) and grey body temperature (right y- 
axis) as a function of redshift. At any redshift slice and for a given 
flux density limit, the red region indicates that all SED shapes of 
the corresponding peak wavelength are detectable, whereas the 
blue and green dashed patterns indicate regions where only 90 
and 70 per cent of the SK07 templates are recovered respectively. 
In the unshaded areas less than 70 per cent of templates are de- 
tectable, reaching zero above a certain redshift. The flux density 
limits used to construct these diagrams are 0.06, 5, 10, 8, 11 and 
13 mJy for the 24, 100, 160, 250, 350 and 500 /.tm bands respec- 
tively, corresponding to the COSMOS 3 cr flux density limits. 
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Figure 4. Selection functions for log (Lif^/LQ)=11.5 using the 
PACS 100-1-160 /im criterion (top panel), the PACS 160 + SPIRE 
250 /.tm criterion (middle panel) and the two used in disjunc- 
tion ([100-1-160 /.im] OR [160-f250Mm]). The plot shows SED peak 
wavelength (left y-axis) and grey body temperature (right y-axis) 
as a function of redshift. At any redshift slice and for a given flux 
density limit, the red region indicates that all SED shapes of the 
corresponding peak wavelength are detectable, whereas the blue 
and green dashed patterns indicate regions where only 90 and 70 
per cent of the SK07 templates are recovered respectively. In the 
unshaded areas less than 70 per cent of templates are detectable, 
reaching zero above a certain redshift. The flux density limits 
used to construct these diagrams are 5, 10 and 8 mJy for the 100, 
160 and 250 /xm bands respectively, corresponding to the 3(t flux 
density limits in COSMOS. The black box in each panel outlines 
the redshift range where all SED shapes with peak wavelengths 
50-140 ^m are detectable — note that the redshift range is larger 
in the 3rd panel. 
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Figure 5. Redshift versus for the Herschel sources (black 
points) in COSMOS, GOODS-S, GOODS-N. Tlie regions below 
the green and blue dashed lines mark the complete parameter 
space derived from the [100+160 ^m] and [160+250 fim] selec- 
tion functions respectively. The red curve corresponds to the 
parameter space traced out by the two criteria in disjunction 
([100+160 Atm] OR [160+250 /^m]), with the hatched red pattern 
indicating the region of incompleteness. In the work presented 
in this paper we use only sources which lie below the red curve, 
which includes any source with 50 < Apeak (m™) < 140, or, equiv- 
alently, temperature within 18 < T(K) < 52. 
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Figure 6. The redshift distribution of the final Herschel sample 
of 1159 sources used in this work; dashed blue line for normal IR 
galaxies (NIRGs), solid black line for LIRGs and red dotted line 
for ULIRGs. 



high redshifts. The small area covered by the red pattern im- 
plies that it is only over a small redshift range that all SED 
shapes are recoverable. However, the large blue and green 
shaded regions indicate that cold SEDs can be detected up 
to high redshift, as long as they do not have a high far-IR 
to mid-IR ratio. 

Ignoring the 24 /^m selection for the moment, Fig. [3] 
shows that our far-IR selection criteria of a detection at 
[100+160 ^m] OR [160+250 /im] wiU resuh in the largest 
number of sources detected in an unbiased part of L-z space. 
This is more clear in Fig. U which shows the [100+160/im] 
selection function (top panel) , the [160+250 fj,m] selection 
function (middle panel) and the two criteria in disjunc- 
tion ([100+160 ^m] OR [160+250 /im]; lower panel) for 
log (Lir/Lq) = 11.5, at the COSMOS limits. The black boxes 
outline the extent in redshift whereby all SEDs with peak 
between 50-140 ^m will be detected. This translates to a 
temperature range of 18-52 K, using the Wien displacement 

law for a ufi, grey body, T{K) ~ , , where h is the 

Planck constant, c is the speed of light in a vacuum and k is 
the Boltzmann constant and we take the dust emissivity (/3) 
to be 1.5. This choice of peak wavelength/temperature range 
is the best compromise between the range of SED types 
probed and the number of objects studied. As we shall see 
in section [6] increasing that range would not have changed 
the measured average properties of the sample but would 
have significantly reduced the statistics. The black box in 
the lower panel of Fig. 3] shows that using the two criteria 
in disjunction, i.e. [100+160 /^m] OR [160+250 /im], allows 
log (Lir/L0)=11.5 sources to be selected up to much higher 
redshifts, than when these criteria are used separately. 

With the aid of the combined selection functions, we 
now define the complete Lir-z parameter space for each sur- 
vey (GOODS-N & S, COSMOS), shown in Fig.[5l where the 
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curves separate the complete (below curve) and incomplete 
(above curve) part of the parameter space. The space be- 
low the curves indicates the Lm-z range, where any source 
with SED peak wavelength within 50 < Apeak (/^rn) < 140, 
or, equivalently, temperature within 18 < T(K) < 52, is de- 
tectable. As mentioned earlier, the [100-1-160 pm] criterion 
(green curve in Fig. [5} , picks up more sources at low redshift 
but excludes more high redshift sources because its selection 
function quickly turns over to favour warm SEDs provid- 
ing a limited unbiased L-z space at high redshift. On the 
other hand the [160-1-250 fim] criterion (blue curve in Fig. [5]) 
performs poorly at low redshift, because the SPIRE 250 j-im 
which mainly drives the combined selection function, largely 
favours cold sources. At z~l, the combination of 160 and 
250 /im turns over as now these bands are sampling 80 and 
125/im respectively, covering the bulk of IR emission. The 
combined criterion of [100-1-160 /im] OR [160-1-250 ^im] (red 
curve in Fig. O is what we thus use to select the final sam- 
ple used in this work. Note that it is the 160 /im band that 
principally drives the combined selection function, whereas 
the PACS 100 fj,m and SPIRE 250 fim in essence provide an 
additional band in the infrared, vital for our analysis. How- 
ever, as they complement each other well, such that most 
SEDs missed at 100 fim are picked up at 250 fim and vice 
versa, the complete region below the red curve in Fig. [5] in- 
cludes many sources which are in the incompleteness regions 
of both the blue and green curves. In addition, this selection 
criterion ensures that the SED peak is well sampled for most 
sources up to z ~ 2. 

Although the selection outlined above is quite con- 
servative since it is unlikely that all SK07 models with 
50< Apeak (mih) <140 are representative of real objects, it 
nevertheless allows us to perform our analysis within a bias- 
free framework with respect to the PACS and SPIRE sur- 
veys. Hereafter, our study concerns only the 1159 sources 
within the complete parameter space below the red curves 
in Fig. [5] The redshift distribution of the final Herschel sam- 
ple is shown in Fig. [S] 

5.2 The 24 /im selection 

In section [5TT] we assembled a complete sample with respect 
to the Herschel bands, which cover the part of the SED pri- 
marily used in our study (see section [6]). It is now important 
to examine whether the requirement for a 24 /xm detection 
affects the completeness of this sample. 

Fig.[3]shows that although the 24 /im and 160 /im selec- 
tion functions cover approximately the same redshift range, 
some SEDs are systematically missed at 24 /im. We investi- 
gate this further by aiming to answer the following ques- 
tions: what types of SEDs are missed, how common are 
sources with such SED types and how does this affect our re- 
sults. The first question is easier to answer. These SEDs are 
ones with high far-to-mid-IR ratio, resulting from a combi- 
nation of parameters in the SK07 formulation, such as high 
extinction (j4v>70) and/or low luminosity and/or low dust 
density within the hot spot, with a detection rate that is 
also redshift dependent. Examples of these SEDs (in the ob- 
served frame) are presented in Fig. [71 and Fig. [8] shows the 
/500//16O - /160//24 colour space these cover; in both fig- 
ures we assume the COSMOS flux limits. Fig. [7] shows that 
such templates are characterised by steep mid-IR continua 
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Figure 7. Example SED templates from the SK07 library, shown 
in the observed frame. In the top panel they are normalised to 
a luminosity of Li^=lO^^ Lq and then scaled and redshifted to 
z=0.3. In the lower panel they are normalised to a luminosity 
of LiR=10^'^ Lq and scaled and redshifted to 2=0.9. At these 
luminosities and redshifts, these templates would not be detected 
down to the 60 /ijy COSMOS 24 /tm flux density limit. However 
they are detected in the Herschel bands at the COSMOS 3 cr 100, 
160, 250 fim limits, also shown in both panels. For comparison, 
the SED of Arp220 (black SED; taken from SK07) is also plotted 
normalised at the appropriate luminosity and redshifted. Note 
that at z=0.9, it is only just detected. 

and strong silicate absorption features at 9.7 and IS/iin, a 
result of high extinction in the SK07 formulation. Note that 
these SED types are easily detected in the Herschel bands. 
In both panels of Fig. |8l we also show the SED of Arp220 
(taken from SK07), normalised at the appropriate luminos- 
ity and redshift. Arp220 is one of the most optically thick 
ULIRGs known (e.g. Papadopoulos et al. 2010), so it is in- 
teresting to examine whether it would be detected in our 
sample. The top panel shows that at low redshift it would 
be detected in all bands, however, at z—0.9 it is only just de- 
tected at 24 /im, but easily detected with Herschel. As this 
SED is redshifted further, it will be missed by the 24 /im 
survey, suggesting that optically thick SEDs, with deep sili- 
cate absorption, would not be in our sample at high redshift 
(^>1)- 

Fig. |8] shows the observed /500//16O versus /160//24 
colours for the COSMOS sample as well as for sources in the 
GOODS samples which faU within the COSMOS complete- 
ness region shown in Fig. (5) The grey-shaded region is the 
parameter space occupied by the SED templates which are 
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Figure 8. Observed /500//16O versus /160//24 colour. The grey 
shaded region shows the colours of SK07 templates, detected in 
the Herschel bands within the complete COSMOS L — z param- 
eter space outlined in Fig. \5\ but missed by the COSMOS 24 ^m 
selection (some examples of these templates are shown in Fig.[7|. 
The black dots are the COSMOS sample and red squares are the 
GOODS-S and GOODS-N sources within the COSMOS L - z 
completeness region (Fig. (5)1. The orange asterisk are the colours 
for the SK07 SED of Arp220 redshiftcd to 2 =1.5. 
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detected in the Herschel bands within the complete Lir — z 
parameter space for COSMOS outlined in Fig. O but are 
missed by the 24 ^m selection. Some overlap between the 
samples and shaded region should be expected, as some 
sources could have high /leo / /24 ratios because of high /iso 
rather than a low /24. However, we see very little overlap, 
suggesting that the detected and non-detected SEDs cover 
a significantly different part of parameter space in terms 
of their far-to-mid-IR colour. This implies that, on average, 
SEDs are missed by the 24 fim selection because of a particu- 
lar feature that reduces the amount of 24 fim observed flux, 
such as a silicate feature, rather than an inherently steep 
far-to-mid-IR continuum. This also becomes obvious from 
the values of /leo / /24 ratios that some templates have: such 
high values of /160//24, up to 6 orders of magnitude, can 
only be caused by a strong 9.7 /^m silicate feature; c.f. with 
what is observed for the Arp220 SED redshifted at z =1.5, 
where the 9.7 //m feature falls in the 24 /im band. In terms 
of the /500//160 ratio, there are only a handful of Herschel 
sources with /500//16O >1, whereas a significant fraction of 
the templates in the grey region have such cold colours. This 
is not surprising as we would expect some of the templates 
that are missed at 24 fim to be overall colder and hence have 
higher /500//160; see some examples in Fig. [T] 

To answer the second question 'how common are these 
SED types?', we compare the GOODS (N and S) to the 
COSMOS colours in Fig. [S] As mentioned earlier, the 
GOODS sources shown in Fig. \E\ are within the COSMOS 
completeness region mapped out in Fig. [5] However, for 
GOODS, the 24 /^m flux density limit is twice as deep as 
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Figure 9. The 24 /.im flux density distribution of the 24 /^m pop- 
ulation in COSMOS, GOODS-N and GOODS-S (green filled-in 
histogram) compared to that of the final Herschel sample before 
(black dashed histogram) and after the Herschel completeness 
criteria are applied (red solid histogram; see section [5. ip . 



it is in COSMOS, so in principle these GOODS sources 
could have SEDs with fi6o/f24 >330. We do not find any 
such sources, in fact we note that GOODS objects have 
fi6o/f24 ratios within the range covered by COSMOS, sug- 
gesting that SEDs with high fi6o/f24 (up to z^2) are rare. 

In Fig. [9] we examine the flux density distribution of 
the 24 fim population in the 3 fields under study, in compar- 
ison to the distribution of Herschel sources before and after 
the Herschel completeness criteria are applied (section |5. If) . 
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Figure 10. Plot of redshift versus total infrared luminosity for 
the local (2 <0.1) samples (see section PS-SII : blue crosses: Clements 
et al. (2010), orange triangles: Hwang et al. (2010) and green 
squares: Buat et al. (2010). The region below the black curve 
marks the complete parameter space corresponding to the IRAS 
1 Jy 60 fim selection function, with the hatched red pattern indi- 
cating the region of incompleteness. Only sources below the curve 
are used in our analysis. 



Note the significant offset between the distributions of Her- 
schel sources and that of the 24 population, suggesting 
that Herschel flux densities are intrinsically correlated with 
bright 24 /xm flux densities. This is unlikely to be an arti- 
fact of our selection, as we are using only 'isolated' 24 
objects, so there is no reason why, in principle, a Herschel 
source cannot be associated with a faint 24 /im source. It is 
immediately obvious from Fig.[9]that the fraction of sources 
which would not be part of our sample because of the re- 
quirement of a 24 jj,m detection is very small. Assuming nor- 
mally distributed flux densities and by calculating the mean 
and standard deviation of each distribution, we can com- 
pute na where n is the number of standard deviations (a), 
at the location of the 24 /^m flux density limit. This gives 
a rough indication of the fraction of sources that are likely 
missed due to the 24 pm selection. Before the Herschel com- 
pleteness criteria are applied, we estimate that 0.1, 1.4 and 
0.01 per cent of sources are unaccounted for in GOODS-S, 
COSMOS and GOODS-N respectively. After these criteria 
are applied, the fraction of missing sources in COSMOS goes 
down to 0.2 per cent, whereas for the GOODS fields it is less 
than 0.06 per cent. Although these are rough estimates and 
rest on the assumption of normally distributed flux densi- 
ties, they do indicate that the fraction of sources missed by 
the 24 selection is very small once our final sample is 
assembled in the complete L — z parameter space using the 
Herschel selection functions (section [SH} . 

In light of this analysis, we conclude that (i) the SEDs 
missed by the 24 /xm selection have high far-to-mid-IR ra- 
tios, mainly as a result of high optical depth in the mid-IR 



and deep silicate features, (ii) these SEDs are not a common 
occurrence amongst IR-luminous galaxies up to z ^2 and 
(iii) our final Herschel sample can be assumed to be com- 
plete in terms of the SED types we can detect. Our findings 
are consistent with the study of Magdis et al. (2011) who 
found the Herschel 24 ^m dropouts to constitute a few per 
cent of the IR-luminous galaxy population at high redshift 
and concluded that they must be sources with stronger sili- 
cate absorption features — see also Roseboom et al. (2010); 
Lutz et al. (2011) for discussion of 24 jum selection effects. Fi- 
nally the third question of 'how would the fraction of sources 
missed affect our results' is discussed at the end of section 
01 



5.3 Local sample selection 

In order to compare the Herschel sample to analogous 
sources in the nearby {z <0.1) Universe, we also assemble 
an IRAS-selected local sample of Lir > 10^" L© galaxies by 
combining sources from Clements et al. (2010), Hwang et al. 
(2010) and Buat et al. (2010). For the Clements et al. (2010) 
objects, we retrieve their published single greybody temper- 
atures and total infrared luminosities calculated using IRAS 
60 and 100 /im data as well as SCUBA data. For the other 
two samples (Hwang et al. 2010 and Buat et al. 2010), we 
use their computed total infrared luminosities and calcu- 
late greybody temperatures by fitting a greybody function 
of emissivity /3=1.5 (see section|3]) to the IRAS and AKARI 
fluxes at A ^60 fj.m. In all cases, the SEDs of the local sources 
have some coverage at ^100 fim, either because of AKARI 
or SCUBA data. 

In order to combine these samples for subsequent anal- 
ysis, we select all sources which are in the complete L — z 
region of the IRAS / 60 fim selection function down to a flux 
density of 1 Jy (see method in section [5] and Symeonidis et 
al. 2011a for the IRAS selection function). Fig. [5]shows the 
curve dividing complete and incomplete parts of parameter 
space. We cut the local samples to include only sources in the 
complete parameter space, where any source with /so >lJy 
and SED peak wavelength within 50 < Apcak(Mni) < 140, 
or, equivalently, temperature within 18 < T{K) < 52, is 
detectable. 



6 RESULTS 

6.1 SED characteristics 

Typical SEDs for the Herschel sample (0.1 < z < 2) 
are shown in Fig. 1111 split into the 3 standard luminosity 
classes: normal IR galaxies (NIRGs; 10^°<Lir<10^^), lumi- 
nous IR galaxies (LIRGs; 10^^ <I/ir<10^^) and ultralumi- 
nous IR galaxies (ULIRGs; 10^^<Lir<10"); see section [3] 
for details on the SED fitting. The Herschel bands cover 
the bulk of the far-IR emission for the majority of sources, 
although in some cases the exact position of the SED peak 
might be underestimated or the slope of the mid-IR con- 
tinuum might not be well constrained due to lack of data 
between 24 nm and the first Herschel band used in the fit- 
ting (Aobs=100 or 160 ^m). Nevertheless, in all cases the 
greybody function covers the bulk of the dust emission giv- 
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Figure 11. Typical SEDs for NIRGs, LIRGs and ULIRGs; y-axis rest-frame luminosity in L©, x-axis rest-frame wavelength (/^m). Black 
points are the available photometry from 24 to 500 ^m. The blue curve is the best-fit SK07 model and the green dotted curve is a single 
temperature greybody fit around the photometric peak in v fv ■ 
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Figure 12. The distribution in (log [Ltot/47r_R^] in units of 
L0 kpc-2) for NIRGS (blue dashed line), LIRGs (black line) and 
ULIRGs (red line). The hatched green region indicates the range 
covered by the SK07 library. The arrows indicate the value of IF 
for M82 and Arp220 (taken from SK07). 



ing a good representation of the average temperature of the 
sources. 

To quantitatively describe the global SED shape we use 
the T parameter defined in section [3l Fig. [12] shows the 
distribution in F for the Herschel sample, split into the 
3 luminosity classes and overlaid on the distribution of all 
templates in the SK07 library. Interestingly, the J- distri- 
butions of NIRGs, LIRGs and ULIRGs show large overlap, 
perhaps surprising as one might expect ULIRGs to exhibit 
a noticeable offset to larger values, simply because they 
are more luminous. However this is not the case, indicat- 
ing that many of the Herschel ULIRGs are described by 
cool/extended rather than warm/compact SEDs, in order 
to reach the same radiation strength per unit area as their 
lower luminosity counterparts. Indeed, the majority of ob- 
jects have 8.5< T <10, suggesting that the IR- luminous 
population up to z ~2 is best described by extended rather 
than compact dust emission. Note that the T distribution of 
the Herschel sample covers only a small range of the avail- 
able parameter space. We find that templates with F >11 
are not representative of any object (within the la uncer- 
tainties on and in fact, only about a 1/3 of the num- 
ber of templates in the SK07 library are representative of 
the sample. This provides useful insight on what SED types 
are observationally confirmed in the context of a physically- 
motivated suite of models. 

Fig. [13] illustrates that the observed distribution in T 
translates to broad-peaked SEDs for the Herschel sample, 
particularly evident when comparing to those of well-studied 
local galaxies such as NGC1808, M82, NGC6240 and Arp220 
(their SEDs all taken from SK07). For the Herschel SEDs, 
the slope on either side of the peak is shallow in antithesis 



Figure 13. Average SEDs for the Herschel sample (black curves) 
for NIRGs, LIRGs and ULIRGs, from bottom to top respectively. 
For comparison we also show the SEDs of Arp220, NGC62040, 
M82 and NGC1808 taken from SK07. 

with the SEDs of the more compact starbursts M82 and 
Arp220, which have J^=10 and 11 respectively. 

6.2 Far-IR colours 

Fig. [14] shows the rest- frame versus I/70/iioo 

colours of the best fit SK07 models to the Herschel sample. 
These bands were chosen for two reasons: (i) they probe a 
part of the SED that is well sampled by our data hence sub- 
stantially constraining the SK07 models in that region and 
(ii) they trace the shape of the SED both around the peak 
(I/7o/iioo) and in the Rayleigh- Jeans side of the continuum 
(^ioo/-f'25o)- For comparison we also include the colours of 
the SK07 library as well as those of nearby galaxies and 
modelled SEDs computed using the GRASIL code (Silva et 
al. 1998; 19990. These are: MlOO, M82, M51, NGC 6946, 
Arp220 and NGC 6090, nearby LIRGs and ULIRGs (Vega et 
al. 2008), modelled colours to represent face-on spirals (Sa, 
Sb and Sc), as well as high redshift gamma-ray burst (GRB) 
host galaxies (Michalowski et al. 2008), in essence young, 
compact, star-forming systems of low metallicity. Nearby 
LIRGs and ULIRGs, such as M82 and Arp220, show warm 
colours, whereas M51, MlOO and M6946 are in the cold part 
of colour-colour space. On the other hand, GRB hosts have 
similar values of I/100/L250 but warmer L70/I/100 colours 
than spiral galaxies. Note that the modelled spirals (Sa, Sb, 
Sc) are located along the edge of the available parameter 
space with colours that become colder down the sequence 
from Sa to Sc. Although the colours of the Sc spiral are 
slightly offset from the parameter space that the SK07 tem- 
plates cover, small discrepancies between SED libraries are 

^ http: / /adlibitum. oat. ts.astro.it/silva/grasil/modlib/modlib. html 
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1 

'-70/^100 

Figure 14. Rest- frame Lioo / ^250 versus L70 /iioo colours for the Herschel sample compared to other star-forming galaxy types indicated 
in the legend. The light grey shaded region represents the colours for the SK07 library, with dark grey shading used for templates of 
9 < < 10 and black shading for templates with > 11. The location of face-on spirals (large white circles with black border) shifts 
from top-right to bottom-left with consecutive morphological class, with Sa being in the top-right. The red writing indicates the position 
of the well-known galaxies MlOO, M82, M51, NGC6946, Arp220 and NGC 6090. 



expected, particularly since there are many input parame- 
ters which can contribute to the final SED shape. 

Overall, the SK07 templates extend over a large range 
in colour-colour space, adequately covering the colours of the 
comparison sample. The SK07 colours show a large spread 
below 1/100/^250 and a narrow tail at large values of 
L100/L250 and Lto/I/ioo- This tail, formed by models with 
J- > 11, is not populated by any of the galaxy groups pre- 
sented here, suggesting that such hot SEDs are not typ- 
ical of any type of star-forming galaxy. The darker grey 
shaded region consists of templates with 9 < < 10 and 
is the parameter space that most Herschel sources occupy 
(see Fig. I12|l . These templates have cool far-IR L1Q0/L250 



colours of <10 and a large spread in Lto/Lioo. A value of 
Lto/I/ioo^I.O ±0.3 ties in with an SED peak between 70 and 
100 fim, the range seen in the Herschel sample (see section 
I6.3|l . The hght grey lower-left part of colour-colour space is 
scarcely occupied by the Herschel sample, indicating that 
very cold, cirrus-dominated SEDs with L100/L250 <3 are 
uncommon in Lir > 10^"^ L0 galaxies. This is not a conse- 
quence of either the selection or the survey flux limits; we 
remind the reader that our sample is complete with respect 
to the SED types that can be probed, hence a deeper IR 
survey is not expected to identify IR-luminous galaxies with 
different SEDs to the ones observed here. Colder colours 
(Lioo/-£'250 <3) might perhaps be more common amongst 
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Figure 15. Dust emperature as a function of total infrared luminosity for the Herschel sample (black points). The green filled circles 
are the mean temperatures computed for 11 Lyu bins (see tabic [2} • The green error bars are the standard deviation in the measured 
temperatures of each bin, whereas the black error bars are the error on the mean. The 2 limiting scenarios of the Stefan-Boltzmann law 
are also shown (solid black lines). The dotted horizontal lines outline the L — T parameter space in which we are complete; see section 
[531 



Table 2. The mean temperature, Icr scatter per Ljr bin and 
error on the mean for the Herschel sample — see Fig. 1151 The 
last column is the median redshift of each Ljj^ bin. 
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more quiescently star-forming galaxies with lower infrared 
luminosities (Lir < 10^° L©), but examining such sources is 
beyond the scope of this paper. 

It is interesting to compare the locus of Herschel sources 
with that of the comparison samples. We see significant 
overlap overall, however many Herschel LIRGs and ULIRGs 
are clearly offset from the region traditionally occupied by 
their local counterparts, displaying colder colours consistent 
with more quiescent star-forming galaxies such as MlOO and 
NGC6090. 



6.3 The luminosity - dust temperature relation 

The dust temperature (T) for the Herschel sample, derived 
as described in section [S] is shown in Fig. [15] as a function of 
total infrared luminosity. Note that although the derivation 
of single average dust temperatures represents simplistic as- 
sumptions with respect to dust properties, such as optical 
depth, emissivity, dust geometry and so on, it is currently 
the only consistent way to characterise and compare statis- 
tically large samples of IR-luminous galaxies over a large 
redshift range. 

The mean temperature of the Herschel sample ranges 
from about 28 to 39 K, increasing with iiR and showing an 
average Icr scatter of 5K; see tabled Our results confirm 
that the choice of temperature range (~18-52 K) over which 
our sample was described as unbiased (section [SJ was ade- 
quate, as we find that the minimum and maximum average 
temperatures are offset by about 10 K from 18 and 52 K re- 
spectively. In fact, IR-luminous galaxies with T <25K and 
T >45K constitute ~6 and ~3 per cent of the total popu- 
lation respectively. 

Since the emission from large dust grains in equilib- 
rium, hence the bulk of the IR emission, is well approxi- 
mated by the black body (or grey body) function, we also 
investigate to what extent we can use the Stefan-Boltzmann 
law, L — AeaT* , to interpret the L — T relation, where A is 
the surface area, e is the emissivity and o is the Stefan con- 
stant. A is proportional to B? , which is in turn proportional 
to the dust mass (Mdust), for constant extinction, so one 
can re- write the Stefan-Boltzmann law as L oc A/dust (or 
L oc R^T*). This spawns two limiting scenarios. The first is 



16 M. Symeonidis et al. 



13.0 



12.5 



12.0 







1 1 .5 



cn 
O 



1 .0 



10.5 



10.0 




4IIIIIIIIIIIIIIUI 




ULIRGs 



NIRGs 



0.5 



1.0 1.5 
redshift 



2.0 



2.5 



Figure 16. The luminosity-redshift parameter space of the final Herschel sample used in the analysis (black triangles), split into bins 
of size 0.2 dex in luminosity and 0.2 in redshift. The bins are outlined in blue for NIRGs, green for LIRGs and orange for ULIRGs and 
further delineated in red if used when computing the average SED shape in that bin (see Fig. 1171 for average SEDs) . 



that the emitting area and/or dust mass is constant which 
would result in an L-T relation of the form: L oc T*. The 
second is that the emitting area and/or dust mass is pro- 
portional to the luminosity [L cc E? or L (x Mdust) with the 
temperature remaining constant for all galaxies. The curves 
representing these scenarios are plotted in Fig. 1151 The ob- 
served L — T relation for the Herschel sample is quite flat — 
2 orders of magnitude increase in luminosity results in only 
a 40 per cent increase in temperature — and hence closer 
to the L oc Mdust hmiting scenario. This suggests that the 
L — T relation is mainly shaped by an increase in dust mass 
and/or IR emitting radius and less so by an increase in dust 
heating. In other words, the average dust temperature of 
ULIRGs is much lower than what one would expect if their 
increased luminosity were the only factor shaping the L~T 
relation. This indicates that the dust masses and/or sizes 
of ULIRGs are larger than those of NIRGs, significantly di- 
luting the effect that their increased luminosity has on the 
temperature. 

Fig. ll6l shows the L — z distribution of the sample, split 
into bins of 0.2 dex in luminosity and 0.2 in redshift. For the 
bins additionally outlined in red (24 in total), containing ^ 5 
objects, we compute the average SED in that bin shown in 
Fig. 1171 Fig. [17] is analogous to Fig. 1161 such that each row 
represents a change in redshift interval, and each column a 



change in luminosity interval. The shaded region represents 
the la scatter around the average SEDs, whereas the red 
SED at the end of each row is the average for that row. 
Consistent with what we observe with regard to the L — T 
relation (Fig. I15|) . Fig. 1 171 demonstrates that there is a shift 
in the SED peak from longer to shorter wavelengths with 
increasing infrared luminosity. This is more clear in the last 
column which shows the average SED for each luminosity 
bin: the SED peak (Apoak) shifts from 86 /im in the lowest 
luminosity bin to 65 ^m in the highest luminosity bin. This is 
also seen in Fig. [I3]where the average SED of each luminos- 
ity class is shown, with the mean and standard deviation in 
Apeak being 86±18 fj,m for NIRGs, 75±18 for LIRGs and 
65±17/im for ULIRGs. Note that the la scatter is large, 
partly because the peak is not always well constrained by 
our data, and partly because of the large diversity in SED 
types (see Figs 1 1 1 1 and 1141 for examples). 

Another feature that appears to change with Lir,, is the 
silicate absorption depth, becoming shallower for higher Lir. 
Our data does not probe the depth of the silicate feature, ex- 
cept over a small window in redshift where it coincides with 
the 24 /xm passband. Hence this trend is likely an artifact 
brought about by model degeneracies. In the SK07 formu- 
lation, visual extinction is tied in to the silicate absorption 
depth, the slope of the mid-IR continuum and the SED peak 
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Figure 17. Average SEDs {vLi, vs Arcst) for the Herschel sample. This figure is analogous to Fig. 1161 such that each row represents a 
change in redshift interval, and each column a change in luminosity interval. The central luminosity and redshift of each bin are shown 
on the left side and top of the plot respectively. Average SEDs are shown for L — z bins with 5 or more objects, outlined in red in Fig. 
1161 The boxes are coloured blue for NIRGs, green for LIRGs and orange for ULIRGs. The solid vertical line in the middle of each box is 
at 100 /im, whereas the dashed line denotes the SED peak. The red SEDs at the end of each row are the average SEDs for that row. 



18 M. Symeonidis et al. 




Figure 18. Left panel: Dust temperature vs L^i for the Herschel sample (green circles) and the local (z <0.1) sample (red squares). 
The mean trends for the two samples are also shown (large filled circles) with the large thin error bars representing the Itr scatter in each 
bin and the short thick error bars representing the error on the mean. Right panel: rest-frame colour (Lgo/Lioo) ^ ^ function of total 
infrared luminosity for the Herschel sample (green circles). The large filled circles represent the mean trend with the large thin error 
bars representing the Icr scatter in each bin and the short thick error bars representing the error on the mean. The solid and dashed 
lines are the local luminosity-colour relation and Ic limits from Chapin, Hughes & Aretxaga (2009). 



wavelength. These quantities are significantly degenerate at 
low redshifts, although for high redshift sources (2 >1), the 
photometry probes further into the mid-IR continuum, plac- 
ing additional constraints on the extinction. Nevertheless, 
although the SED shape in the near/mid-IR is more reli- 
ably reproduced for the high redshift sources, the observed 
trend of decreasing silicate depth with increasing Lir is most 
likely artificial. 

6.4 Evolution in dust conditions 

In Fig. 1181 we compare the properties of the Herschel sources 
to those of the local sample, assembled as described in sec- 
tion [S3] The local luminosity-temperature and luminosity- 
colour (C; L(io/L\m) relations, the latter in functional form 
from Chapin, Hughes & Aretxaga (2009), are shown in Fig. 
1181 left and right panels respectively. The Leo/iioo colour 
has been used extensively to characterise the dust tempera- 
ture of local samples and analysis of /7?AS'-selected galaxies 
has shown that more luminous sources have higher colour 
temperatures than their less luminous counterparts (e.g. 
Dunne et al. 2000; Dale et al. 2001; Dale & Helou 2002; 
Chapman et al. 2003). For both L — T and L — C relations 
we note a systematic difference between the Herschel and lo- 
cal samples, with the former displaying lower values of T and 
ieo/iioo- This can be interpreted as evidence for evolution: 
high redshift IR-luminous galaxies have more emission long- 
ward of ~60 pim compared to their low redshift analogues, 
lowering the average dust temperature. This is consistent 
with results from section 16.21 Fig. 1141 where we see that the 
Herschel sample extends to colder far-IR colours than local 
LIRGs and ULIRGs. 

Recent results from the Planck collaboration (Part 16) 
on the dust properties of nearby /7?j45-selected sources 
showed that many local 10^° < Lir < lO^'^ galaxies ex- 
tended to lower temperatures than previously reported. This 
does not affect our work as (i) within the scatter and given 



that emissivity is a free parameter in their fitting, their mea- 
sured dust temperature in the 10^*^ < Lir < 10^^ luminosity 
range are consistent with the ones we report here for the lo- 
cal sample (Fig. I18p and (ii) as described in section 15.31 we 
assemble the local sample in an unbiased part of parameter 
space. Furthermore, lack of additional submm/mm data for 
our nearby sources, would not change the average dust tem- 
peratures measured, as these trace the peak dust emission 
and are thus insensitive to inclusion of photometry signifi- 
cantly longward of the peak (see also Magnelli et al. 2012). 

Although a systematic reduction in the dust tempera- 
ture and colour of IR-luminous galaxies from the local to 
the high redshift Universe is observed, it is worth examin- 
ing how this translates to changes in the SED shape with 
redshift within the Herschel sample. In Fig. 1171 we see some 
evidence for a shift in the SED peak to longer wavelengths 
along the rows, i.e. with increasing redshift, for example in 
the log LiR=11.4-11.6 bin. However, overall, it is not clear 
how the SED shape evolves with redshift. This is likely a con- 
sequence of sparse sampling of the SK07 templates by our 
photometry, exacerbated by the small redshift range covered 
along each row, statistical uncertainties in each bin as well 
as the fact that in many cases, sources do not cover each 
L — z bin uniformly (see Fig. I16|) . A different way to inves- 
tigate evolution in dust properties is to repeat this exercise 
with our computed dust temperatures, in order to remove 
model-dependent uncertainties, although the other sources 
of uncertainty outlined above would remain. Fig. ll9l shows a 
2-D image of the dust temperature of IR luminous galaxies 
(local and Herschel samples combined) as a function of red- 
shift and luminosity, in bins of 0.2 and 0.2 dex respectively. 
Besides an increase in average dust temperature vertically 
along the luminosity axis, consistent with our analysis in 
section 16.31 on the L — T relation, we also note an overall 
reduction in the mean dust temperature, horizontally, along 
the redshift axis. This is more pronounced when considering 
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Figure 19. 2-D image of tlie L — T — z space for IR-Iuminous galaxies from the local Universe to z = 2. This includes both the local 
and the Herschel samples. Objects are divided into Ljjj^ (y-axis) and z bins (x-axis), and the average temperature of each bin is shown 
as a greyscale intensity map. The colourbar on the right is the temperature key for the map. The temperature bins extend from 22.5 to 
45 K. White colour indicates unpopulated or underpopulated (i.e. <5 objects) L ~ T — z bins. 
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Figure 20. The fraction of T<35K NIRGs (orange squares), LIRGs (black diamonds) and ULIRGs (red triangles) as a function of 
redshift. The local sample is at z ^O.! and the remaining bins include sources from the Herschel sample only. 



the first and last bins of each row, however in some cases it 
is also evident along the length of the row. 

Considering the above trends, it is interesting to ex- 
amine whether the evolution in dust temperature we ob- 
serve, is luminosity dependent, i.e. whether dust conditions 
of ULIRGs evolve at a different rate to those of NIRGs 
(Fig. [20)l . We place the separation between cold and warm 
at T=35K, corresponding to the mean temperature for the 
Herschel sample shown in Fig. 1181 Fig. I20l shows an increase 
in the fraction of cold ULIRGs with redshift, from 5 per cent 



in the local Universe to about 30 per cent at z ~l-2. Also 
the fraction of cold LIRGs increases from about 60 per cent 
locally to about 80 per cent at high redshift. However, as for 
Figs [17] and 1191 in some bins, the L — z parameter space is 
not sampled uniformly. The higher redshift bins in Fig. [201 
include a larger fraction of more luminous and hence warmer 
sources, implying that the fraction of cold sources is likely 
underestimated, resulting in the observed downturn of the 
computed fraction. 

Note that given the depth of our data, currently this 
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is the best attainable coverage of tiie L — z plane in the 
horizontal (z) direction within the unbiased framework we 
define in this work. This does not mean that LIRGs at z > 1 
and ULIRGs at z > 2 will not be detected by Herschel or 
other facilities, rather it implies that it is not currently pos- 
sible to measure the aggregate properties of the LIRG and 
ULIRG population at those redshifts. On the other hand, 
with larger area Herschel surveys we expect to cover the 
gap between the local and high redshift Universe in the ver- 
tical (L) direction, enabling better sampling of the L — z 
plane, where large survey area is needed, and hence achieve 
better statistics for logLia/L© >11.5 galaxies up to z = 1. 

With respect to the selection at 24 /im, we saw that we 
are likely missing a few per cent of SEDs with steep far- 
to-mid-IR continua and deep silicate features, about half of 
which have /500//16O >1 extending to higher values than we 
find in the Herschel sample (see Fig. [8] and section 15. 2p . 
This suggests that many of the sources missed by the 24 /im 
criterion might have lower temperatures than the average 
temperatures derived for the Herschel sample. Although this 
fraction of sources is very small and hence unlikely to change 
our results, it would nevertheless only serve to strengthen 
the differences we observe between the local and high red- 
shift sample. 



7 DISCUSSION 

7.1 The properties of the IR-luminous population 

The Herschel sample under study consists mainly of LIRGs 
(64 per cent), with ULIRGs constituting 20 per cent, con- 
sistent with what is expected in the redshift range probed 
(0.1 < z < 2), whereas the remaining 16 per cent are nor- 
mal IR galaxies (NIRGs; 10^° < Lm < lO" Lq). The dust 
temperatures of the sample show large scatter from 15 to 
55 K, however the mean temperature ranges from 28±4 to 
39±6 K, with ULIRGs being on average about 10 K warmer 
than NIRGs. Similarly, we see a shift in the average SED 
peak wavelength from 86±18/im for NIRGs to 65±17/im 
for ULIRGs. The sample is best described by cool/extended, 
rather than warm/compact SEDs, and broad peaks, trans- 
lating to low values of IF, a parameter which we defined as a 
measure of the overall IR SED shape. In addition, there is a 
large overlap in colour-colour (1/100/^250 ~ L70/L100) space 
between the Herschel sample and other star-forming galaxy 
types, the coldest of which (spirals and young compact 
star- forming galaxies) are at I/100/I/250 ^7 and the warmest 
(ULIRGs and starbursts) at I/ioo/-i'250 ^7. For the Her- 
schel sample, we noted a roughly equal number of sources 
above and below Lioo/I/250 ^^7. It is worth mentioning that 
only about a 1/3 of the SK07 templates are representative 
of the Herschel sample. The Lioo/i25o>20, Lioo/i'25o<3 
and I/7o/I/ioo>2 regions of the SK07 library are scarsely 
populated. Moreover, the distribution in of the SK07 li- 
brary extends to much higher values (J-"~15), than what 
is observed in the Herschel sample (J^<11). This indicates 
that compact, hot-dust-dominated SEDs or very cold cirrus- 
dominated SEDs are not typical of the IR-luminous popu- 
lation. 

How do our findings compare to other studies of IR 
galaxies? As mentioned earlier, we have performed a rig- 
orous analysis of selection effects in order to minimise any 



biases which would interfere with our results. Hence, due 
to the nature of our study, we are sensitive to most, if not 
all, IR galaxy types with Lir > 10^" L0 at z=0.1-2. Conse- 
quently, results on the properties of IR galaxies from previ- 
ous studies with Spitzer (e.g. Symeonidis et al. 2009; Kar- 
taltepe et al. 2010a; Patel et al. 2011) and SCUBA (e.g. Ko- 
vacs et al. 2006; Coppin et al. 2008; Santini et al. 2010) are 
all within the parameter space we probe (see also Magnelli 
et al. 2012). The results reported recently using Herschel 
data (e.g. Rowan-Robinson et al. 2010; Hwang et al. 2010; 
Smith et al. 2012) are also within the parameter space we de- 
fine here for the IR-luminous population. However, we note 
that although cold, cirrus-dominated SEDs such as those 
reported in Rowan-Robinson et al. (2010) and Smith et al. 
(2012) are part of our sample, they represent a small frac- 
tion (<6 per cent) of the IR-luminous population and are 
not the prevalent SED types. 



7.2 The L - T relation 

The increase in dust temperature as a function of infrared 
luminosity we observe here, is similar to the L — T relation 
that has emerged from most infrared population studies (e.g. 
Dunne et al. 2000; Dale et al. 2001; Dale & Helou 2002; 
Chapman et al. 2003). Since the emission from large dust 
grains in equilibrium, hence the bulk of the IR emission, is 
well approximated by a black (or grey) body function, we 
aimed to understand the L—T relation within the framework 
of the Stefan- Boltzmann law. As described in section [631 the 
two limiting scenarios of the Stefan-Boltzmann law are: (i) 
L oc r* with R (the radius of the emitting region), or Mdust 
(the dust mass) kept constant and (ii) L oc B? or L oc Mdust 
with T kept constant. The former scenario would produce a 
steep L — T relation, whereas for the latter the L ~T rela- 
tion would be flat, with the increase in luminosity tying in 
with an increase in surface area of the emitting body or the 
dust mass. We flnd that the L — T relation for the Herschel 
sample lies closer to the latter scenario suggesting that its 
shape is mainly driven by an increase in dust mass or extent 
of dust emitting region and less so by energetics. In other 
words, it seems that the increased dust heating in ULIRGs 
is diluted by an increase in their physical size and/or dust 
mass, such that the L — T relation is diverted away from 
the L oc r* scenario. However, although more strongly star- 
forming galaxies such as ULIRGs are predominantly more 
massive (e.g. Dave 2008; Shapley 2011) and dustier systems 
(e.g. Magdis et al. 2012) with warmer average dust temper- 
atures, we find that their SED shapes are not substantially 
different to their lower luminosity counterparts. In particu- 
lar, we do not see a clear segregation in far-IR colour-colour 
space or in the range covered by as a function of infrared 
luminosity. Interestingly, this suggests that properties such 
as optical depth, dust extinction, extent of IR emitting re- 
gions etc., which determine the overall SED shape are not 
substantially different between high redshift ULIRGs and 
their lower luminosity systems. 



7.3 Evolution of the IR-luminous population 

In this work we found evidence that the dust temperatures of 
Herschel sources are systematically colder than equivalently 
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luminous galaxies in the local Universe. The rigorous analy- 
sis of survey biases we performed ensures the validity of this 
result, as our final sample selection was sensitive to all IR 
galaxies with 18 < T(K) < 52. We found sources at z >0.5 
and log [Lir/Lq] >11.0 to be on average between 5 and 10 K 
colder than their z <0.1 counterparts. We also noted a sys- 
tematic offset to colder rest frame Leo /Lioo , Lioo /L250 and 
Lro/iioo colours for the Herschel sample in comparison to 
local equivalent sources. We believe that this is unlikely to 
be due to an increase in extinction, which by removing flux 
from the mid-IR and adding it in the far-IR could mimic a 
lower dust temperature. As discussed earlier, sources with 
high extinction and deep silicate absorption are missed by 
our 24 ^m selection at z >1. However, they do not consti- 
tute more than a few per cent of the population (see also 
Magdis et al. 2010). Moreover, for the high redshift sources, 
where the Herschel photometry can more accurately con- 
strain the slope of the mid-IR-to-far-IR continuum, we find 
that the sample is mainly fit with low extinction (shallow 
silicate absorption depth) models. 

Note that although previous studies have shown that IR 
galaxies colder than local equivalents exist in abundance at 
high redshift — e.g. results from 750 (e.g. Rowan- Robinson 
et al. 2005), SCUBA (Kovacs et al. 2006; Pope et al. 2006; 
Coppin et al. 2008), Spitzer (e.g. Symeonidis et al. 2008; 
2009; Kartahepe et al. 2010a; Patel et al. 2011), BLAST 
(Muzzin et al. 2010), Herschel (e.g. Hwang et al. 2010, 
Rowan-Robinson et al. 2010; Smith et al. 2012) — for the 
first time we determine that the mean dust temperature 
of the IR-luminous population as a whole decreases as a 
function of redshift. Moreover, we find that the decrease in 
temperature is also a function of infrared luminosity, i.e. 
the temperatures of more luminous objects show a stronger 
decline from the local to the early Universe. We note that 
almost all NIRGs, up to z ~0.5, have temperatures below 
35 K, whereas for LIRGs the local cold (T <35K) fraction 
is 60 per cent increasing to about 90 per cent at z '-^0.5. 
The ULIRGs show the largest increase in the cold galaxy 
fraction, from about 5 per cent at z < 0.1 to 30 per cent 
at 2 = 1 — 2. This is interesting as it implies that LIRGs 
undergo more modest evolution than ULIRGs, the former 
showing a 2-fold increase in the fraction of cold galaxies, 
whereas the latter a 6-fold. Moreover, the cold LIRG frac- 
tion in the local Universe is about 12 times higher than the 
cold ULIRG fraction, however, we see that this difference 
decreases to about 2.5 at z ~0.8. Similarly the cold NIRG 
fraction in the local Universe is about 50 per cent higher 
than the cold LIRG fraction, whereas at z ^GA these frac- 
tions are about the same. This is evidence that cold galaxies 
become more dominant at high redshift. However, it also in- 
dicates that there might be a lower limit in the average dust 
temperature of the IR-luminous population, at T ~ 25 K, to- 
wards which systems tend. This is not to say that T <25 K 
IR-luminous galaxies do not exist, but these would be at the 
tail of the temperature distribution. 

As discussed earlier, the L — T relation would be com- 
pletely flat, e.g. all IR luminous galaxies would have an av- 
erage temperature of ~25 K, if their sizes or dust masses 
increased in proportion to their total IR luminosities (see 
also Fig. [T^ . As it stands, this is not the case, and so 
the increased luminosity succeeds in heating up the dust 
to a higher temperature. However, the decrease in average 



dust temperature with redshift suggests that high redshift 
LIRGs/ULIRGs must have more extended IR emitting re- 
gions and/or higher dust masses relative to their lower red- 
shift counterparts, causing the L — T relation at high redshift 
to become flatter than the local one. Described phenomeno- 
logically, we observe that the temperature evolution of IR- 
luminous galaxies is more rapid if their local temperatures 
are much higher than 25 K, such as for ULIRGs, than if they 
are closer to 25 K such as for NIRGs. Our results are in agree- 
ment with the work presented in Dunne et al. (2011) who 
flnd strong evolution in the dust mass density, proposing 
that IR-luminous galaxies are dustier at z ~0.5 compared 
to today, corresponding to a factor 4-5 increase in the dust 
masses of the most massive galaxies. Moreover, CO mea- 
surements support the idea of extended instead of compact 
star-formation in high redshift ULIRGs, which have >>kpc 
CO sizes, in contrast to local equivalent sources which are 
more concentrated (<lkpc) (e.g. Tacconi et al. 2006, lono 
et al. 2009; Rujopakarn et al. 2011). Moreover, signiflcantly 
higher gas fractions in z ~ 1 disc galaxies than in nearby 
discs have been reported (e.g. Tacconi et al. 2010). Our re- 
sults also agree with Seymour et al. (2010), who studied 
the comoving IR luminosity density (IRLD) as a function 
of temperature, proposing that cold galaxies dominate the 
IRLD across < z < 1 and are thus likely to be the main 
driver behind the increase in SFR density up to z '^l. 

Although the IR-luminous population and particularly 
LIRGs and ULIRGs seem to have more extended dust dis- 
tribution and/or higher dust masses at high redshift, the 
analysis presented here cannot constrain whether these sys- 
tems are characterised by a merger-induced or isolated star- 
formation history or where they are located in the star- 
formation rate-stellar mass parameter space (e.g. Reddy et 
al. 2006, 2010; Wuyts et al. 2011; Whitaker et al. 2012; Zahid 
et al. 2013). Morphological studies of IR-luminous galaxies 
have presented contrasting results — e.g. Bell et al. 2005; 
Lotz et al. 2006 report that more than half of LIRGs at z> 
0.7 are gas-rich isolated spirals, whereas other studies (e.g. 
Le Fevre et al. 2000, Cassata et al. 2005, Bridge et al. 2007; 
de Ravel et al. 2009) claim an increase in the merger rate out 
to z ~1 suggesting that more than half of IR luminosity den- 
sity out to z ~1 is a result of some merger event. Moreover, 
Zamojski et al. (2011) and Kartaltepe et al. (2012) flnd more 
than 70 per cent of ULIRGs up to z ^2 to be mergers. These 
findings are hard to reconcile and recent evidence from Lotz 
et al. (2011) shows that these differences might simply be 
the result of increasing gas fractions at high redshift, as the 
timescale for observing a galaxy to be asymmetric increases 
in tandem with the gas fraction, with the resulting dust ob- 
scuration also being a key factor. With the work presented 
here, we are not able to test this argument nor examine the 
morphological evolution (if any) of IR galaxies. Neverthe- 
less, our description of the IR-luminous population is inde- 
pendent of morphological classification. Our results indicate 
that the gas-rich environment in the early Universe might 
have set or enabled different initial conditions in these sys- 
tems, resulting in the observed differences between IR galax- 
ies at high redshift and their local counterparts. The increase 
in cold (T <35K) galaxy fraction reported here suggests a 
greater diversity in the IR population at high redshift, par- 
ticularly for (U)LIRGs. In contreist, the dust properties of 
the local (U)LIRG population are more uniform and as a 
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result they are not archetypal of the (U)LIRG population 
as a whole. 



8 SUMMARY & CONCLUSIONS 

We have examined the dust properties and infrared SEDs of 
a sample of 1159 infrared-selected galaxies at ^=0.1-2, using 
data from Herschel /P ACS and SPIRE and Spitzer/MIPS 
(24 ^m) in the COSMOS and GOODS fields. The unique 
angle of this work has been the rigorous analysis of survey 
selection effects enabling us work within a framework al- 
most entirely free of selection biases. The results we thus re- 
port should be considered as representative of the aggregate 
properties of the star-formation-dominated, IR-luminous 
(LiR>10^'^) population up to z ^2. 
We conclude that: 

• IR-lumiuous galaxies have mean dust temperatures be- 
tween 25 and 45 K, with T < 25 and T > 45K sources 
being rare. They are characterised by broad-peaked and 
cool/extended, rather than warm/compact SEDs, however 
very cold cirrus-dominated SEDs are rare occurrences, with 
most sources having SED types between those of warm star- 
bursts such as M82 and cool spirals such as M51. 

• The IR luminous population follows a luminosity- 
temperature (L — T) relation, where the more luminous 
sources have up to 10 K higher dust temperatures. However, 
the effect of increased dust heating is not solely responsible 
for shape of the L — T relation. We find that the increase 
in dust mass and /or extent of dust distribution of the more 
luminous sources dilutes the increased dust heating, flatten- 
ing the L — T relation and driving it towards the limiting 
scenario of L oc J?^ or L oc Mdust with T=const. 

• High redshift IR-luminous galaxies are on average 
colder with a temperature difference that increases as a func- 
tion of total infrared luminosity and reaches a maximum of 
10 K. For the more luminous (log [Lir/Lo]>11.5) sources, 
the L — T relation is flatter at high redshift than in the 
local Universe, suggesting an increase in the sizes and/or 
dust-masses of these systems compared to their local coun- 
terparts. 

• The fraction of T<35 K gala:xies increases with redshift 
as a function of total infrared luminosity. Although NIRG 

fractions arc consistent in the local and high redshift Uni- 
verse, LIRGs show a 2-fold increase and ULIRGs a 6-fold in- 
crease. This suggests a greater diversity in the IR-luminous 
population at high redshift, particularly for ULIRGs. 
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APPENDIX A: REDSHIFTS 

Here we examine the reliability of the photometric redshift 
catalogues used in this work (see section 12. 2p for sources 
which satisfy the initial criterion for the selection of our 
sample (section [2TT}: 24 ^m sources that have detections (at 
least 3cr) at [100 and 160 /xm] OR [160 and 250 ^m]. 

For GOODS-N, we use the redshift catalogue of Berta et 
al. (2011) which includes spectroscopic redshifts assembled 
from various sources and photometric redshifts compiled us- 
ing the EAzY code (Brammer et al. 2008) and up to 14 pho- 
tometric bands. Berta et al. (2011) note that the fraction of 
outliers, defined as objects having Zspcc-Zphot/(l+2spcc)>0.2, 
is about 2 per cent for sources with a PACS detection. For 
more details we refer the reader to Berta et al. (2011). In 
Fig. lAll (top-panel) we compare the photometric and spec- 
troscopic redshifts for the GOODS-N sample, indicating the 
10 and 20 per cent uncertainty regions. There is excellent 
agreement in the redshifts. 

For GOODS-S, we use the MUSYC redshift catalogue 
of Cardamone et al. (2010), complemented with the MU- 
SIC redshift catalogue of Santini et al. (2009). Most red- 
shifts come from Cardamone et al. (2010) with the addi- 
tion of some from Santini et al. (2009). The Cardamone 
et al. (2010) photometric redshifts are compiled using up 
to 28 photometric bands and the EAzY code (Brammer et 
al. 2008), resulting to high accuracy Zspoc-Zphot/(l + 2spec) 
~0. 00822 out to high redshift. For more details we refer 
the reader to the Cardamone et al. (2010) paper. Fig. lAll 
(middle panel) shows a comparison of photometric to spec- 
troscopic redshifts for the sources in the GOODS-S sample 
for both Santini et al. (2009) and Cardamone et al. (2010) 
catalogues. There is excellent agreement between the spec- 
troscopic and photometric redshift in both catalogues, with 
the catastrophic failures, which we define as a more than 20 
per cent offset in (l-f^), being <1 per cent of the spectro- 
scopic sample. 

For COSMOS the spectroscopic redshifts are from Lilly 
et al. (2009) and the photometric redshifts from Ilbert et al. 
(2008), derived using up to 30 photometric bands. The es- 
timated accuracy of a median Zspec-Zphot/(l + 2spec)=0.007 
for the galaxies brighter than i = 22.5. For the sources which 
host X-ray detected AGN we substitute the Ilbert et al. red- 
shift with a photometric redshift from Salvato et al. (2009); 
(2011). These are derived with a combination of AGN and 
galaxy templates and are hence more accurate for galaxies 
hosting AGN. Fig. lAll (lower panel) shows a comparison 
of photometric to spectroscopic redshifts for the COSMOS 
sample sources. The agreement is again excellent, with only 
about 1 per cent of sources outside the 20 per cent uncer- 
tainty region. 

We next examine the redshift distribution of our final 
Herschel sample used in this work assembled in section (5] 
This is shown in lA2l for the 3 fields. Note that the photomet- 
ric and spectroscopic redshift distributions agree, although 
the photometric redshift distribution tails off to higher red- 
shifts. We examine whether the use of photometric redshifts 
would have any effect on our results, by reproducing one 
of our main figures using only spectroscopic redshifts (Fig. 
IA3|I . We see that using only spectroscopic redshifts does not 
change the overall differences we find between the Herschel 



and local sample, however it does significantly reduce the 
statistics. 
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Figure Al. Comparison of spectroscopic and photometric red- 
shifts for sources in the GOODS-N, GOODS-S and COSMOS 
samples which satisfy the initial criterion for the selection of our 
sample (section 12.11 1 : 24 lira sources that have detections (at least 
So-) at [100 and 160 /im] OR [160 and 250 /im]. The dashed lines 
are the 10 per cent boundaries and the dotted lines are the 20 
per cent boundaries. In the middle panel, the asterisks are the 
Cardamone et al. (2010) photometric redshifts and the diamonds 
the Santini et al. (2009) photometric redshifts. 
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Figure A2. The redshift distribution of the final Herschel sample 
used in this work in COSMOS, GOODS-S and GOODS-N. The 
filled in histogram shows the total redshift distribution whereas 
the red histogram represents the spectroscopic redshifts and the 
blue dashed histogram the photometric redshifts. 
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Figure A3. This figure is reproduced from the main part of 
the paper, but solely with spectroscopic redshifts. We see that 
the overall results remain the same, however the statistics are 
reduced. 



